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A truncated moment formalism for general relativistic radiation hydrodynamics, based 
on the Thome's moment formalism, is derived. The fluid rest frame is chosen to be the 
fiducial frame for defining the radiation moments. Then, zeroth-, first-, and second-rank 
radiation moments are defined from the distribution function with a physically reasonable 
assumption for it in the optically thin and thick limits. The source terms are written, focusing 
specifically on the neutrino transfer and neglecting higher harmonic angular dependence of 
the reaction angle. Finally, basic equations for a truncated moment formalism for general 
relativistic radiation hydrodynamics in a closed covariant form are derived assuming a closure 
relation among the radiation stress tensor, energy density, and energy flux, and a variable 
Eddington factor, which works well. 

§1. Introduction 

Radiation fields and their interaction with matter often play a crucial role in 
many astrophysical contexts. For example, the critical roles of photon pressure dur- 
ing proto-star and massive-star formation and of neutrino heating and cooling in 
supernova core collapse and explosion are well-known among many other phenom- 
ena. To theoretically clarify these phenomena, it is necessary to solve hydrodynamic 
equations as well as radiation transfer equations. For strictly handling the radiation 
transfer, it is necessary to numerically solve the Boltzmann equation, taking into 
account the absorption, emission, and scattering terms. However, this equation has 
3-|-3-|-l dimensional form (3 dimensions in real and phase spaces, respectively, and 
1 dimension in time), and furthermore, the time scale for the interaction between 
matter and radiation is often shorter than the dynamical time scale of the system. 
Thus, in the current computational resources, it is not feasible to perform a well- 
resolved numerical simulation with a sufficient grid resolution. A certain approximate 
method incorporating key features of radiation effects is often required in numerical 
astrophysics. In particular, no useful formalism for multi-dimensional simulation in 
general relativity has been well developed (but see Refs. [T])-[3])). Note that in spheri- 
cal symmetry, this equation is simplified to a l-|-2-|-l dimensional form as formulated 
in Ref. Simulations with similar formalism were performed in Ref. [5]), and sub- 
sequently, sophisticated simulations including the state-of-the-art microphysics were 
achieved, e.g., in Refs. ED ,[7]). However, the effort has been paid only to the spherical 
symmetric simulations. In this paper, we derive an approximate formalism of radi- 
ation hydrodynamics in general relativity, in which a numerical simulation will be 
feasible capturing the physically important ingredients. 



typeset using PTPTeX.cIs (Ver.0.9) 



2 



Historically, a popular method for approximate radiation hydrodynamics is a 
flux-limited diffusion (FLD) methodP In this method, the radiation flux density is 
in general assumed to be described by the radiation energy density, and resulting 
evolution equation for the radiation energy density becomes a diffusion-type equa- 
tion in the optically thick region (cf. § 5). In this case, the propagation speed of 
characteristics may be larger than the speed of light, although in general relativ- 
ity, the causality must not be violated. Another drawback of the FLD scheme is 
associated with the presence of constraint equations (Hamiltonian and momentum 
constraints ) in the initial value problem of general relativity: In numerical relativity 
for multi-dimensional problems, we usually solve the evolution equations of Einstein's 
equation and matter equations self-consistently. As a result, the constraints are sat- 
isfied within a numerical error. However, in the case that we do not solve the energy 
and momentum equations for the radiation field self-consistently, the constraints are 
violated. In the FLD method, one solves an equation only for the radiation energy 
density component, and hence, the constraints will be violated in general. 

Truncated moment formalisms have been also proposed for an approximate solu- 
tion of radiation hydro dynamics .'^''^ In this approach, one derives a set of covariant 
equations for multi-pole moments defined from the distribution function of radiation. 
Then, assuming that higher-order moments may be neglected and imposing closure 
relations, a closed covariant form of basic equations is derived. With an appropriate 
choice of the closure relation, the causal relation can be preserved, and furthermore, a 
solution of the radiation transfer in the optically thick and thin limits can be derived 
from the resulting equations. In this paper, we derive a truncated moment formal- 
ism in general relativity following the covariant formalism developed by Thome's 
In addition, we derive a closed coordinate-independent formalism including the ab- 
sorption, emission, and collision terms, focusing specifically on neutrino transfer in 
high-density and high-temperature medium. 

The paper is organized as follows: In § 2, we review the covariant moment 
formalism derived by Thorne.'^l' In § 3, a truncated moment formalism is presented, 
assuming a physically reasonable specific form for the distribution function. In § 4, 
source terms of the moment formalism are written in terms only of the radiation 
field variables employed in our truncated moment formalism, focusing specifically 
on neutrino transfer. In § 5, approximate solutions for the radiation fields in the 
optically thick limit are derived. In § 6, we propose a closure relation among the 
radiation scalar, vector, and tensor. We also derive the characteristic propagation 
speeds of the radiation field in the optically thick and thin limits for the chosen 
closure relation. In § 7, hydrodynamic equations coupled with the radiation fields 
are derived. In § 8, radiation hydrodynamic equations in a slow-motion limit (usually 
referred to as Newtonian radiation hydrodynamic equations) are derived. Section 9 
is devoted to a summary. Throughout this paper, Greek (a, /3, 7 • • • ) and Latin (i, 
j, k - ■ ■) subscripts denote the spacetime and space components, except for which 
always denotes the angular frequency of radiation (which never be the subscript of 
space or time), always denotes spacetime coordinates. We assume to use the 
Cartesian coordinates as the spatial coordinates for simplicity. Unless otherwise 
stated, the units of c = 1 = /i are used, where c is the speed of light and h the Planck 
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constant, denotes the Boltzmann constant. 

§2. Moment formalism of Thorne 

First, we review the Thome's moment formahsm.-' In the first step, he defines 
an unprojected moment of massless particles associated with a moving medium as 

M(°/"--°'=(x/^) = I lMl2^f^^^y^p'^.p'^2 . ..p'-.dv;, (2-1) 

where / is the distribution function of the relevant radiation, z^' = —u^p'^ the fre- 
quency of the radiation in the rest-frame of the medium (i.e, in the rest-frame of 
the fiducial observer) with being medium's four velocity, p^^ the four-momentum 
of the radiation, and dVp the invariant integration element on the light cone, k, 
here, is positive integer, 1, 2, • • • . As pointed out by Thorne,'^ the choice of the 
fiducial observer is crucial when deriving a good truncated formalism from his mo- 
ment formalism. In the following, the fiuid, coupled with the radiation, is chosen 
as the medium .13 'flSJ Namely, the frequency, v, in M^"^^"^'"'^* always denote the 
frequency measured in the rest-frame of the fluid throughout this paper. This choice 
is crucially helpful when computing the source terms of the radiation equations. 

We note that it is possible to choose any fiducial frame in the moment formalism. 
However, we have to keep in mind that for a truncated moment formalism in a closed 
form, it is necessary to assume a closure relation which is determined by a physically 
reasonable assumption. In the dense medium, radiation is strongly coupled to the 
matter field. This implies that at the zeroth order, the radiation is in equilibrium 
with the medium, and radiation fiow (measured by an observer comoving with the 
matter) is a small correction. To reproduce this feature in the closure relation, the 
best method seems to choose the fiuid rest frame as the fiducial frame. 

We also note the following: As a result of our choice of the fiducial frame, the 
argument frequency in the distribution function is always the frequency measured in 
the fluid rest frame. By contrast, the argument variable should be in general the 
frequency in the laboratory frame (although any frame can be taken), if one fully 
solves the Boltzmann equation that the distribution function obeys. 

The Boltzmann equation is written in the form^J 

where S denotes a source term and r the affine parameter of a trajectory of radiation 
particles. In any orthonormal frame, the invariant integration element is given b>0 

where is the four-momentum of the radiation in the local orthonormal frame. In 
the local rest frame of an observer comoving with the fluid. 



dVp = vdvdQ, 



(24) 
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where J dQ denotes integrations over solid angle on an unit sphere. 
For the following, we write in the form 

i^=p" = ^(u" + r), (2-5) 

ar 

where is a unit normal four- vector orthogonal to u"; iai.^ = 1 and = 0. 

Using this decomposition of p", Eq. ()2-ip is rewritten to give 

^a^a2-ak = „^ j /(z., X^^) (u^l + ) (u^^ + T^' (2-6) 

Here, the angular dependence is included in 1°' and the following relations hold. 



dnf 



0= f dm^i^p, — f dQei^ = 

J ^7tJ 3 

^ j dnn^e^i^ = (^h'^^h^^ + h'^^h^^ + h'^^h'^^^ . (2-7) 

hafj is the projection operator defined by 

hal3 ■■= gal3 + U^Ui^. (2-8) 

Following Thorne,'^ we denote jVf^°^i°2 M^^'' . Taking the covariant derivatives 
of M^^''^ , we obtain a covariant equation with respect to real-space coordinated 
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where denotes the covariant derivative associated with the spacetime metric ga/B , 
and 



^^^"=1,^ J 5(z^,i7,x'^,/)(n"i +^«2)...(^afe (2-10) 

Here, the spacetime derivative is taken holding and the frequency derivative is 
taken holding spacetime location. It should be noted that Eq. ()2-9p has a coordinate- 
independent form as stressed by Thorne.'^' Also, the following relation is worthy to 
note: 

^(.tV = -^(.t- (2-11) 

Thus, the rank-(fc + 1) equations include the lower-rank equations. 

Since the frequency, u, in Eq. ()2-9p denotes the frequency observed in a fluid- 



rest frame, not in the laboratory frame, M^^'^ is not directly related to the spectrum 
observed in the laboratory frame. However, if the fluid is assumed to be at rest 
in a distant zone far away from a radiation source where we observe the spectrum, 
the radiation moments in the fluid-rest frame agree with those in the laboratory 
frame. We suppose that the present formalism will be used for the system that 
this assumption holds, e.g., supernova stellar core collapse and merger of compact 
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objects. Thus, it is possible to directly compute the radiation spectrum from M^^'', 
if we estimate it for r — )• cxd (cf. Appendix A for an example). 

Integrating Eq. ()2-9p by z^, we obtain (for each species of the radiation compo- 
nent) 

V^M^*^ -{k- l)M^''^^V^up = (2-12) 

where 

di^M^^^" and S^" = J duS^^^K (2-13) 

Equation (|2-12p is essentially the same as the moment formalism derived by Anderson 
and Spiegel.'^ We note that the second-rank tensor M"^ is equal to the energy- 
momentum tensor for one of the radiation components. 

In the following, we analyze only the second-rank part of Eq. ()2-9p . truncating 
the higher-rank parts (in § 5, we partly use the third-rank equation for deriving a 
solution in the absence of closure relation). In the next section, we develop such 
formalism. 

§3. Truncated moment formalism 

First of all, we define the following moments: 

J^,y.= iy^ J f{u,[2,x^')d[2, (3-1) 
H(^y.= v^ j rf{u,Q,xndQ, (3-2) 
L^^^ := 1/3 j n^fiu, Q, x'')dQ, (3-3) 
AT^^/^T := J.3 j ('^iPfff^jy^ xt^yQ_ (3.4) 

Here, all these integrals are assumed to be performed in the local rest frame comoving 
with the fluid, and v denotes the angular frequency of radiation measured in this 
local rest frame. The second- and third-rank moments are denoted by 

^(.T = -^H^"^^ + ^(.>^ + ^(.>" + ^(.t' (3-5) 

+ ^(?-" + h'^y + ^('^" + ^(^r- (3-6) 

and the total stress-energy tensor for the radiation is 

/•oo 

where the summation denotes to sum up for all the species of the radiation. 
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In our truncated formalism, (i) we formally define the zeroth-, first-, second- and 
third-rank moments from the distribution function, and (ii) we solve the evolution 
equations only for the zeroth- and first-rank moments. For the optically thick region, 
this is approximately equivalent to assuming that the degree of anisotropy of the 
distribution function in the fluid local rest frame is weak and that the distribution 
function is approximated by 



f{u, xn = fo{u, + ffiu, X^ia + ff{y, X^iaip. (3-8) 

Here, /o, /f , and /^'^ do not depend on the propagation angle of radiation in the fluid 
local rest frame and /^^ is a traceless tensor with respect to hajs (i.e., f2^hap = 0). 
We assume that |/o| is much larger than the absolute magnitude of /f and /^^. For 
the expansion of Eq. (|3-8p , we obtain 



J(,) = ATiv'h, (3-9) 
47r 



H,:. = — z^Vf , (3-10) 



= f -^(/o^"^ + \ff) = + ^-Vf , (3-11) 

^(T = \ {Hi>P' + + ' (3-12) 



where we used the relations (|2-7p . Thus, /o and /f are directly related to J(i^) and 
H^^y and /^^ to the traceless part of Lf^^ , respectively. Because of the truncated 

expansion for f(y, [2,x^), we naturally obtain a closure relation for N^^'^^^ . 

For the optically thin limit, by contrast, we should first give a physical assump- 
tion in the laboratory frame because the radiation does not interact with matter. 
We employ the assumptions that the radiation should propagate with the speed of 
light and the radiation flow at each spacetime point should be pointed to a null 
direction. The former assumption then implies that the radiation flow is pointed 
to a null direction in any frame (although the spacetime coordinate basis changes). 
Thus, for such region, the distribution function may be written in the form (see ^6.11 
for details) 

f{u, f2, x^") = 47r/f (i/, x'')6{[2 - Qi), (3-13) 

where f2f denotes the flow direction in the fluid rest frame, and /f (z^, x'^) is the partial 
distribution function of = i7f . Then, the radiation moments are calculated to give 

J(,) = W/f, (3-14) 
F^^) = ATTu'hef, (3-15) 

L(^f =47rz.Vf^f^f, (3-16) 

N^;^f^ = A7ru^fieie^£], (3-17) 

where if denotes the unit vector of the flow direction (observed in the fluid-rest 
frame). In Appendix A, we illustrate that the assumption of (|3-14|) - (|3-17|) would 
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be appropriate for providing the radiation field solution in the optically thin-limit 
medium. 

The equations for J(^) and H^^^ are derived from the second-rank part of Eq. 



as 



where 



The acceleration a" and the shear are defined by 

a" := u'^Vpu'^, (3-20) 



VyUs + VsUy . (3-21) 



To obtain a closed set of the equations, we have to determine L^'^!^ [N^^^^^"^ is given by 
Eq. ()3T2p or ()3T7|) ]. Instead of solving the equation for this, which may be derived 
from the moment equation of third rank, we will assume a closure relation for it; an 
artificial (but physically reasonable) relation between L^'^!^ and {J (jj-^ , H ^^■^) will be 
assumed (see ^|6]). 

Substituting Eq. p-5p into Eq. p-lSp . the evolution equations for J(,^) and H^^^ 
are obtained as 



V/3Q(,7 + Q(:;)V^n"-^[KL(,7n^ + iV(,7^)V;3n,]J =%.5(,"), (3-23) 



V„Q(^) + Q^2';Vpu^ - V^n„)] = (3-22) 

/i^^ "/^ _L n Y7 _ A [„r r "T,,/? _L at "/^t 

where 

Q(^) := -M(^f = J(,)^.- + H^^^, (3-24) 

■= ^"^^(Jf = ^(")^' + ^(?- (3-25) 
The frequency-integrated equations are 

V^Q" + Q"^V;3U„ = (3-26) 
/ifca(V/3Q"^ + Q'^V^sn") = (3-27) 

where 

/"OO /"CO 

Q-:=y^ dvQ^^y Q''^-=J^ d^Q^:,)- (3-28) 

Thus, the equations are not in the conservation form; the reason is that and Q^^ 
are not conservative quantities even in the absence of the source terms. 
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Instead of using Eq. (j3-5|) . M^^^ may be written by 

M^^l^ = E^^^n'^n^ + F^^n^ + F^^^n^ + (3-29) 

where n" is a unit vector orthogonal to the spacehke hypersurface. -^(j^; ^^^d 

P^^^ may be regarded as radiation fields measured in the laboratory frame. We note 
again that the meaning of the frequency, v, is unchanged; it is the frequency observed 
in the fluid rest frame. To obtain the quantities fully defined in the laboratory frame, 
we need the transformation of v to the frequency measured in the laboratory frame. 
However, in the moment formalism, we do not consider such transformation, as 
already mentioned. 

We however assume that = in the far region with r — )• oo as mentioned 
in §2. Our primary purpose is to develop an approximate formalism which can be 
used for simulation of stellar core collapse and merger of binary compact objects. 
For such purpose, this assumption is acceptable. Thus, for r — )• oo, we suppose 
that Ei^y-^ = J(^^-^, Ffy^ = H^'j^y and Pf^^-^ = ^(^) ■> ^"^^ ^ agrees with the frequency 
measured in the laboratory frame. Therefore, if Ey is extracted in a distant zone in 
the numerical simulation, we can obtain the spectrum of the radiation. 

In the 3+1 formulation of general relativity, 

n"=(-,-^), (3-30) 

where a is the lapse function and the shift vector. Then, E(^^y ^{u)^ ^iy) 
are defined by 

= Mh'^«^/3, = n„7,^ Piij,=M^i:i^, (3-31) 

where 7q,^ is the three metric 

7a/3 := 9al3 + naUji- (3-32) 

Because Ffy^ria = Pf^^-^Ua = 0, we have the relations Ff^^-^ = Pf^^^ = 0. 

Here, we consider a formalism in which E'^j^) and P^^-^ are evolved, and P^y^ is 
determined by a closure relation. Then, J^^) and -ff^"^ are determined by 

J{u) = E{u)W^ - '^F{y)WUk + P^^-^UiUj, (3-33) 
= {E^^,)W - F^^^UkWpn^ + wh^pF^i;^ - h%u,P^^y (3-34) 

where w = au^ . We note h°'pn^ = n° — wu°' and riah^^Jisk = —wuk- For the later 
convenience, we give relations 

Q{u)'^a = -J{u)W + Hf^^^ria = -E(^^)W + F^^^Uk, (3-35) 
<3(")7ai = J{v)Ui + = wF^i - P(^^fiUk. (3-36) 
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As mentioned before, it is natural to assume that Uj = (to = 1/q ~ 1) for the dis- 
tant zone far from the radiation source. Then, both the asymptotic power spectrum 
densities, £'(^) and J(^), agree with each other, because the frequency v agrees with 
that in the laboratory frame. 

The evolution equations for Ef^^-^ and are written in the conservative forms 

as 



d 



d 



where 7 is the determinant of 7jj and Kij the extrinsic curvature. 
The frequency-integrated equations are 

= a^[P'^Kij - F^dj Ina - 5"n„], 



^\ - Edta + FkdiP'' + ^P^'^dajk + aS^li, 



where 



E :-- 



duE, 



dvFi 



dvp;{. 



(3-37) 



(3-38) 



(3-39) 
(3-40) 

(3-41) 



Equations ()3-39p and ()3-40p have fully conservative forms, because E and Fi are the 
conservative quantities in the absence of the source terms and gravitational fields. 
Thus in the numerical simulation, it will be better to adopt basic equations based 
on Eqs. (|3^ and (f3^ . 



§4. Source terms 



The source terms for the second-rank radiation field equations, S^^y have to be 

written in terms of the radiation moments {J(^^^, H L f^^) . 5^"^ is formally written 

as 

S(^^ = j (I?, x^) (n" + r (4-1) 

where B(^^^ is the so-called collision integral. In the following, we assume that S(^^^ 
and B(^y-^ are written as a function of the phase-space coordinate defined in the local 
rest frame of the fluid. The real coordinate, x^, is arbitrarily chosen. 

We derive the source terms focusing on the neutrino transfer in a high-density 
and high-temperature medium. We show that under a reasonable and often-used 
assumption (anisotropy of the collision integral is small), the source terms are totally 
written in terms of J^y), H,^., and -^("f • 



10 



4.1. Neutrino absorption and emission 

First, we consider the absorption and emission of neutrinos by nucleons and 
heavy nuclei such as n+Vf, o p+e~ , p+D^ <-)• n+e+, and (Z, A)+i'e -H- {Z—1, A)+e~ , 
where n, p, e^, Ve {^e), and {Z,A) denote neutrons, protons, electrons (positrons), 
electron neutrinos (anti neutrinos), and heavy nuclei, respectively. For these cases, 
the collision integral is written in the forrrPJ'HZ) 

= - /(^, ^?,^^)] - (4-2) 

where denotes the emissivity, A(,y) is the neutrino absorption mean free path, and 
/(z^, n,x^) denotes the distribution function of relevant neutrinos (in the following, 
we omit the argument in /). j(^) and A(,^) are quantities independent of the 
neutrino propagation angle, f2. 

The integral of Eq. (|4-1|) is easily performed to give 

5(.")= 47rj(,)Z.3„- - ( J(,)..- + (j(,) + A(-;p 

= (JM + [(Jiu) - Jm)^'' - H^u)] , (4-3) 

where 

J5> := 4vrz.3_JH ^ 4vr^V^^(^), (4-4) 
and f^'^{i') is the equilibrium distribution function. For neutrinos (fermions), 

where Hc and T are the chemical potential and the temperature for the corresponding 
species of neutrinos which are in thermal equilibrium with matter. For the following, 
we define the opacity as 



K{u) ■■= + >^(^J), (4-6) 



and thus. 



Q a 



(J(^,^)-J(.))^" . (4-7) 



4.2. Neutrino- electron scattering 

Neutrinos are scattered by electrons, nucleons, and heavy nuclei. In the case of 
electron scattering, the collision integral is generally written a^^''^ 

S(^) = j v'^dv'dQ'[f{v', Q'){1 - f{v, Q)]R'''{v, v\ w) 

-/(z., n){\ - fii^', n')}R°-\i^, u', u;)], (4-8) 

where oj is the cosine of the scattering angle, and R}^ and are the scattering 
kernels. Following Refs. [TT|). I12|) . we approximate these kernels by taking the terms 
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up to the linear order in to, i.e., 



(4-9) 
(4-10) 



uj is related to the angular part of the momentum-space coordinates Q = {0,f) and 
f2' = {9', if') of the ingoing and outgoing neutrinos by 



bj = COS 6* COS 9' + sin sin 0' cos (99 — if'), 
and thus, we have the following relations 

Consequently, the collision integral is written as 



(4-11) 



(4-12) 



(4-13) 



and the source term is 
dv' 



Q a. 



{(47ri/3 - J, 



+ 



("' 



where -^("j^ is the traceless part of L^^"!^: 



a/3 

(y) 



(4-14) 



(4-15) 



Because of the approximation in which the terms up to the linear order in lo for the 
scattering kernel is taken in account, the source term is totally written in terms of 
J(u), ^("), and L^'^^. 

4.3. Pair production 

For the thermal neutrino pair production and pair annihilation, the collision 
integral has the following forrcP^-IUll 

S(,)= J u''du'dn'[{l- f{u,f2)}{l- f{iy',n')}R'P'^°{u,u',u) 

- f{u, n)f{u', Q')R-^'\v, v\ uj)] , (4-16) 
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where / and / are the distribution functions of neutrinos and anti-neutrinos, re- 
spectively, and i^P''" and R^^'^ are the integration kernels for pair production and 
annihilation, respectively. This integral can be performed in the same manner as in 
the neutrino-electron scattering: For the expansion up to the linear order in w, 

RP'"{u, u', uj) = Rl'\u, u') + Rf°{v, v')uj, (4-17) 
i?^'^°(i^,zy',L<j) = Rl''''{u,u') + Rl"'^{u,v')uj, (4-18) 



we obtain 



dv' 



- ( J(,)^z- + J(,,)iig-(^, I.')] , (4-19) 

where the quantities with bar denote the radiation moments for anti-neutrinos. 
Again, the source term is totally written in terms of J(u)j ^{°)' ^^"^ ' 

4.4. Isoenergy neutrino scattering 

In the neutrino scattering with nucleons and heavy nuclei, the energy exchange 
may be assumed to be zero. In such isoenergetic neutrino scattering, the collision 
integral is written a^nj-O 



[ dn'[f{i^,n') - f{iy,n)]R''°{iy,uj), (4-2o) 



where v denotes the angular frequency of the ingoing and outgoing neutrinos. Fol- 
lowing Refs. [TT]). ll2p . we again approximate the kernel R^^°{i',uj) by taking the terms 
up to the linear order in lo as 



Then, 



and thus, 



where 



B, 



Q OL 



ISO TT a 



R^i^) - ^RTi^) 



(4-21) 
(4-22) 
(4-23) 
(4-24) 



Therefore, the source term depends only on the radiation flux (not on J(,y)), as clearly 
shown in Ref. [I2D [cf. Eq. (A.46) of Ref- IT^ ]. 
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§5. Optically thick limit 

In this section, we derive a solution of the radiation moment equations for the 
hmit that the radiation is optically thick to absorption and scattering with matter 
fields. For neutrinos, such a limit may be realized only for a special phenomena, 
such as stellar core collapse and merger of binary neutron stars, in which an extreme 
state of high density and high temperature is likely to be realized. However, it is 
quite useful to derive the solution for the idealized situation and to confirm that the 
derived solution is physical, for validating a new formalism. 

In the following, we analyze the radiation moment equations derived in the 
previous sections, taking into account only the neutrino absorption, emission, and 
isoenergy scattering for simplicity, for which the source terms of the collision integral 
are written by Eqs. ()4-3p and ()4-23p . Then, the source terms of Eqs. (|3-22p and (|3-23p 
are in the form 

hkaSf^^-j = (5-2) 

where 

= «H + 'tf°)- (5-3) 
We define a mean free path as the inverse of the opacity 

which measures the optical thickness of neutrinos: In the optically thick limit, /(^,) 
is much smaller than a characteristic length scale of system (e.g., a stellar radius 
or p/lii^V fj,p\), and the radiation fields may be expanded by /(j,) assuming that it is 
sufficiently smalP (this is the so-called Thomas approximatiorP'Ef). Assuming that 
is also of order Z(j^), we can expand the radiation moments as 

JM = Jlu)+lMJiu)+0{ll)), (5-5) 
H^^^ = l^^^H^^^- + 0{ll^), (5-6) 

^(.t = ^^M^"' + V)^!;r + 0(^(.))' (5-7) 

where /(j^) = . In the following, we determine J(j^) , i?^ , and L^^^ up to the first 

order in /(^) (^{u))- ^or this purpose (specifically for deriving it is necessary to 

analyze not only the second-rank moment equations but also the third-rank moment 
equations 

V,M(^f ^ - M^^f^'Vsu, - ^{i^M-f^'Vsu,) = S-^^. (5-9) 
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Thus, we have to take into account the fourth-rank moment as 

jja^f^^S ,_ ^3 J l»lPp(Sj^^^ (5.10) 

and expand up to the zeroth order in Z(^) as 

^ a/375 ^ J_ jeq (;ja/3/^57 + /j^^/j/^T + /j«7/i/35) + 0(/(^)). (5-11) 

In addition, we have to calculate the second-rank moment of the source term which 
is 



Q a/3 



1 



(5-12) 



where /^"^-j = = + 47rz/^_Rg°. 

The expanded solutions for the radiation moments are determined from the 
expanded equations for Eqs. ()3-22p and (|3-23p in each order in If^^y Their zeroth- 
order equations give the first-order solutions 



J, 



(1) 



1 r 



il)a 



3 L 
3 



d 



(5-13) 
(5-14) 



and -^^l^j"^ is derived from the zeroth-order part of the third-rank moment equation. 
Taking into account that this term is traceless and its component is perpendicular 
to n°, we obtain 



where 



L 



151 



4j;q _ ^{vJ':\) 

{") du ^'^r 



a 



al3 



(5-15) 



a 



(5-16) 



For the frequency- integrated case, all these solutions agree with those in Ref. [T]). 
The physical meaning of the first-order correction for the frequency-dependent equa- 
tions (except for the diffusion effect associated with the term ^ aJ^^^) is essentially 
the same as that for the frequency-integrated case: The first-order corrections of 
J(^), H^^y and L^^ are associated with the fluid expansion {O = VqU'^), the fluid 
acceleration {ap = li^VaU^), and the fluid shear {aap)- 

We note that from the derived first-order solutions, the second-order solutions 
for J(j^) and i/^"^ are easily derived from Eqs. (|3-22p and (j3-23p : 



J, 



(2) 

(^) 



r(l)" 



15 



d 



dul2, (^) 



Llf \ CLfy -\ i-^l \ (Jn 



(5-17) 



H, 



(2)a 



d_ 
dv 



where J^^^j and H^f^^ are the coefficients of /^^^ and f^^^ , respectively. For clarifying 
the order and for simplicity, we here assume that l(^y^ and /(^) are constants. 

It is interesting to note that for the frequency-integrated caseP the magnitude 
of the first-order terms is proportional to J'^°' where for each species of neutrinos 



r(2)" 



\7i _i_ 4 A'^'>„'y 



)}], (5-18) 



jeq 



(5-19) 



For the frequency-dependent equations, the magnitude depends universally on 



ld_ 



du 



hv 



(5-20) 



Because of the presence of a factor 1 — f^'^{i') and f'^°'{v), these first-order corrections 
play a role only by neutrinos of energy around /Uc — ky^T ^ hu <^ + l^hT, i.e., near 
the Fermi surface, if the corresponding species of neutrinos is degenerate. This is a 
reasonable consequence and characteristic property for fermions. O 

The first-order solutions for J(^), H^^y and L^"^ may be used to constitute a 
diffusion equation from which the first-order solutions are produced. Substituting 
the first-order solution into Eq. (|3-22p with replacement of J"^^-^ to J(i,) gives 



u d [ 



where 



J, 



J: 



15 



(5-21) 



(5-22) 



Thus in general, the diffusion equation is modified by the acceleration and shear 
motion of the fluid and by neutrinos, if we do not assume the slow motion of the 
fluid. 



*' We note that only electron neutrinos can be degenerate in general. For anti electron neutrinos, 
fic < and for muon and tan neutrinos, = 0, when these neutrinos are in thermal equilibrium 
with matter. Thus, these are not degenerate in general. 
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We note that in an FLD approximation, the first-order solution for H^^-^ is mod- 
ified as 



If 



H,''. = ^ ^4 

and is then substituted in Eq. ()3-22p . With this prescription, the equation for J(,^) 

reduces to a wave equation with the characteristic speed ~ c for the case that /(j,) is 
much longer than a characteristic length scale of the system. 

§6. Closure relations 

In the truncated moment formalism derived in §3, we proposed to solve the 
equations for Ef^^,^ and F^^^^ but not to solve that for P^J^y which is assumed to 

be determined in terms of E'^^) and F^^y In this section, we propose a physically 
reasonable closure relation. 

6.1. Optically thin case 

In the limit that the optical depth is zero, the emission, absorption, and scat- 
tering are negligible. When the source term of the radiation field equations can be 
neglected, the radiation freely propagates, and the radiation moments should obey 
a wave equation with no source. 

One example for such region is the asymptotically fiat region, far from the ra- 
diation source where curved spacetime effects as well as hydrodynamic effects play 
a tiny role (e.g., we may consider ~ n'^, J(i^) ~ -^(i^)' ^{^) ^ ^{t^) ^ ^-li'sady 
mentioned in §2 and 3). Thus, any closure relation assumed has to satisfy at least 
the equations in the fiat spacetime. 

For the fiat spacetime, we obtain the equation for F^J^ from Eq. (j3-37|) 

5.(^/^^h) = 0, (6-1) 

where r] is the determinant of the fiat three metric r]ij . This provides a reasonable 
solution of -F^J^ for the spatial infinity; for the spherically symmetric fiow, F^^^ oc r~^, 
and for the plane symmetric fiow, F^^-^ =constant for the flow direction. On the other 
hand, Eq. ([3^ gives 

dkiVvPiu) j) = ^Pil)djmk. (6-2) 

For an appropriate solution of F(^), the following closure relation is the first 
candidate (and is that we finally choose): 

This choice satisfies Eq. (j6-2p in the asymptotically flat region. This is regarded 
as a general relativistic extension of the so-called Ml closur^^^'fi^li for the optically 
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thin-limit region. In this case, E(^,y-j oc r ^ for the spherical symmetric flow and 
=constant for the plane symmetric flow, and hence, the reasonable condition 

is guaranteed. Furthermore, P^^^^jk is equal to E(^^-) in this condition, guaranteeing 

the necessary condition for the radiation fields, ga/BT^^^ = 0. 

It should be noted that this choice with no modification may not be accepted in 
general relativity, because two of the characteristic speeds of the wave equations for 
and F^^^ may exceed the speed of light (see § 16.41 and Appendix B). The pure 
choice of this closure relation is allowed only for 



^H = V^^.^H^H' (6-4) 

in which the characteristic speed is guaranteed to be equal to the speed of light. This 

is guaranteed in the optically thin limit. However, for > ^ ^ijF^^-^F^^-^ which 

may often happen in a not-completely free streaming region, the characteristic speed 
may exceed the speed of light. This implies that an appropriate modification in the 
grey region is required in the choice of this closure relation to satisfy Eq. ()6-4p (see 
§ 16.31 for a candidate choice of variable Eddington factor in the grey region and 
Appendix C for a satisfactory test result). 
Another possible candidate is 

(^y^_ (6-5) 



This choice also satisfies Eq. (|6-2|) in the asymptotic region. With Eq. (|6-5|) . Eq. 
for the spherical and plane-symmetric stationary flows is written as 

VvF^I^^Okfi' = 0, (6-6) 
where is a unit vector parallel to F^^^.^ 

For the spherical and plane-symmetric stationary flows, the condition (|6-6|) is guar- 
anteed, and hence, the closure relation (|6-5p is acceptable. 

With this choice, the characteristic speed of the equation for is approx- 
imately the speed of light in the asymptotically flat region (see § 16. 4|) . However, 
P^i^ljk is not a priori guaranteed to be equal to Ei^^y, for this condition to be satis- 



fled, y jijF^^^F^^^ has to be equal to E'f^^) but the condition will not be satisfled in 
the near zone, i.e., near the emission source (cf. Appendix B). Only in the far zone, 
this condition seems to be followed from Eq. (j3-37p because the radiation propagates 
with the speed of light as \F^^-^\ ~ -£-(1/)- Because of this reason, we choose Eq. ()6-3p 
for the closure relation. 
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It is reasonable to suppose that Eqs. ()6-3p and (|6-4p are satisfied in the optically 
thin hmit. Because we have F^-^ = E(^y^f^ where /" is a unit spatial vector, = 

1, and orthogonal to n", n'^fa = 0. J(,^), H^^y and L^^ are rewritten by 

= M^^^^UaU/B = E^^^){u'^qaf, (6-8) 
Hf^:^ = -M^^^^uph^^ = -E(,)U^qph^^q\ (6-9) 
L(^"f = Mj^'h^^hf = E^^^h^^hfq^', (6-10) 

where q" = n" + is a null vector. Defining 4:Ttv^ = E^y'^{u"'qa)'^ and = 
—hp^q^/{u^q^), we obtain Eqs. ()3-14p - p-16p . and find it reasonable to assume 
Eq. (|3-13p as the distribution function in the optically thin limit. It is easy to 
confirm that IfUa = 1 for the above definition. 

Finally, we have to give a closure relation for the third-rank moment. We propose 
to employ 

iV,f^= '^^^^'^(-"''(-f(-\ (6-11) 



or 



M-f^ = ^±J^±J^. (6-12) 

a/s {v) {u) 

Here is related to ^{1)^ and P^^^ by Eq. (|3-34p . Associated with the choice 
of Eq. (li^D, we choose Eq. (lfrTT|) . 



6.2. Optically thick case 

As shown in § O in the optically thick limit with l^^jj-^ — )• (or K(^) — )• cxd), -Z^^^j' 



IS written as 



1,.«. 4 



(6-13) 



The correction term of 0(/(j,)) is associated with the so-called radiation viscosity. 
With this prescription, we can incorporate the first-order effect of L^'^^ in /(j^) without 

solving the third-rank moment equation. (Note that N^'^^^ should be given by 
Eq. (|3-12p in the present formalism.) For the case that velocity of the medium is 
much smaller than the speed of light, \v^\ <^ c, the term of 0(/(j^)) may be neglected. 
This term will play a role for the medium moving around a black hole, such as in 
black hole accretion disks. 

For numerical computation, we have to transform the relation of Eq. ()6-13p to 
the relation of P^^j as a function of and F^^^y For this purpose, we omit the 

term, d{uJ{^y'^)/dv^ for simplicity. The reason is that in its presence, is not 
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written by and F^J^^ in a straightforward manner (although it is possible to do 
in an approximate manner). For the frequency-integrated case, this term vanishes, 
and hence, we may say that the radiation viscosity effect is taken into account in 
a frequency- averaged way. However, in this treatment, the low-energy neutrinos, 
which should not contribute to the radiation viscosity for degenerate neutrinos, may 
incorrectly play a role. To avoid this unphysical contribution, it will be appropriate to 
artificially reduce /(j^) to zero for hu fic — k\^T when treating degenerate neutrinos. 

Assuming that Eq. (j6-13p holds with the omission of the third term, we have 
the relations 

r4u;2 -1-1 
E{u) = ^ o-Q J(u) + 2-f^(i/)j^-', 



Hi 



L3 



wui + ai 



Ui 



where = 'y^^Uj {Vi = Ui), and 



4/ 



(Jo = -rr-o- ' naUji, ai 
15 



AI 



15 



(6-14) 
(6-15) 

(6-16) 



Also, we used H^^^^^u"' = and H^^^ = {aw) ^ H^^j^-^iV' . Equations (|6T4p and (|6-15p 
constitute simultaneous equations for J(j^) and H(^,y^i. Inverting them yields 

-2w^ 



J, 



1 1-1 

— + O-Q 



{2w'' - 1)^(,) - 2wF^^^Uk 



(6-17) 



H -If 1 



+ [{Aw'^ + l)ui + 6wai + 3aoUi]Ff^l^^Uk 



. (6-18) 



Note that F^^^Uk = F^^,)^^ but H^^.-^kV^ + H^^u^; H^^^ = (7'='-/3S'-^z^/aw;)F(,);. 



Also wao = —aiV. Then, P/-? is given by 



p ij 



J, 



3 15 ^^^^ 



+ H^l^V^+H^l^V\ (6-19) 



where J(^) and i^(^)(= 7'V^(^)) are given by Eqs. (|6T7)l and (l6T8)l . With this 

closure relation for P(y^, the necessary condition for the radiation fields, QapT^^ = 0, 
is guaranteed to be satisfied. 

We note that with the closure relation ()6-19p . the first-order term in /(,^) may be 
accidentally larger than the zeroth-order term for a high value of fXy. Thus, it may 
be necessary to change the definition of /(j^) as 



mm 



1 ^ ( y'^'^k \V2 



(6-20) 



where is a coefficient smaller than unity. 
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6.3. Grey zone 

For a solution of the radiation fields in the optically grey zone, in general, it is 
necessary to fully solve the radiation transfer equation in general relativity. However, 
it is not possible in the framework of truncated moment formalism and far beyond 
the scope of this paper. We propose an approximate method which is essentially the 
same as the variable Eddington factor method.^^' In this prescription, P^*^^ is given 

by 

^(^f = ^(^'(.pthin + ^^^(^(Jpthick, (6-21) 

where x is the so-called variable Eddington factor, which is x = 1/3 in the optically 
thick limit and x = 1 in the optically thin limit. Following Ref. [T3|l . we choose that 
X is a function of F, for which in general relativity, the candidates are 



7„-F,\F/' \ 1/2 



F:=[^^^f^] , (6-22) 



and 



F:=(^^^^^) . (6-23) 

For the optically thick and thin limits, F = ^ and F = 1, respectively. For giving a 
correct value of F in the optically thick limit, Eq. ()6-23p should be chosen because 
i?^"^ should be zero in the comoving frame; if the fluid has a large uniform velocity, 
the value of F in Eq. (j6-22p would be highly different from zero even in an optically 
thick medium. For giving a correct value of F in the optically thin limit, both 
Eqs. ()6-22p and (j6-23p can be chosen, because in such a limit, M^^^ is proportional 

to J[u)P'^P^ (p" is a null vector) and F = 1 for the null fluid in both definitions (see 
§3). For this reason, we choose Eq. (j6-23p for F. 

With the choice of (|6-23p , F obeys an algebraic equation for a given set of E(^^^ 
and F,\. This can be written in the form 



F^ = (6-24) 



where for , Eq. ()3-29p is used with Eq. (|6-2ip . In numerical simulation, we have 
to solve this equation numerically. 

Livermore proposed several functions for x{P)^ e.g., 

3 + 4F2 

X = , 6-25 

5 + 2V4 - 3F2 ^ ' 

In the Appendix C, we employ Eq. ()6-25p . and show that it is likely to work well. 
However, it should be kept in mind that it might not be the best one and better 
closure relations should be further explored. 
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For completeness, we have to provide N^^^^^ . We propose to employ 



afi'y 



3x-l 



thin + 



3(1 - X) a/37 



(^) 



/thick) 



(6-26) 



where (A''^ ")''''^)thin and (A'"^ ")^''')thick are given by Eqs. ()6-lip and ()3-12p . respectively. 



6.4. Characteristic speed 

For numerical computation with conservation schemes, it is necessary to know 
characteristic speeds. Furthermore, the analysis of the characteristic speed is helpful 
to check whether the proposed closure relation is acceptable or not (i.e., it is smaller 
than the speed of light). 

The characteristic speed of the radiation moment equations is computed from the 
Jacobian matrix (e.g., Refs. [TC|) - [T7|) ). For the conservative variables E^^^-^ and F(^ly 
the Jacobian matrix for the x direction is (in the following, we omit the subscript v) 



A. 



ab 



a- 



a 



a- 



dE 
dP% 

'dE 



07 



07 
a- 



dP 



dP- 

a 



y 



a 



a- 



dPx 
dP% 



dP% 
-/3^ + a- ^ 



a- 



dP% 



dp, 



a- 



dPz 
dP% 



dP% 



(6-27) 



The characteristic speeds are computed from 

det{Aab - Xlab) = 0, 

where lab is the unit matrix. 

For the optically thin case with the closure relation 

px 



A 



± a- 



and with the closure relation (16-5 



and A = -jS"" + aE 



PuPk 



(double) , 



A = -/3^' and A = -/3^' + a- 



(triple). 



(6-28) 



(6-29) 



(6-30) 



As we already pointed out, |A| can exceed the speed of light for the closure relation 
(j6-3p if the opacity is not zero limit; \EF^ / F}^P^\ may exceed the speed of light. 
Hence, it is not allowed to be employed without appropriate choice of the variable 
Eddington factor for the optically grey zone. By contrast, with the closure relation 
(j6-5p . the characteristic speed is smaller than the speed of light (but in this case, the 
tracefree condition for the stress-energy tensor is not satisfied in general, as already 
mentioned) . 

For the optically thick limit with = 0, 



A = -/3^ + 
and A = — 



2w^p ± ^a2^xx(-2^t;2 + 1) - 2wV 
2u;2 + 1 
(double). 



(6-31) 
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where p = aV^/w. For w = 1 (p = 0), the first one reduces to 



X = -l3-±a^^—, (6-32) 

and thus, we obtain a well-known characteristic speed for the radiation fluid in the 
diffusion limit (~ l/\/3). For w ^ oo (p — )■ a/ sj^xx < 1); 

A^-/3^'+p, (6-33) 

and thus, A approaches to a local light speed (but never exceeds it). For /(^.^ ^ 0, 
the characteristic speed is modified in a complicated form. However as far as /(^,) is 
small, this effect is not important. 

The general formula for the closure relation (|6-3p is written in the following 
manner. For i-direction, 

X = -/3'±a ,^ and \ = -P' + aE-^ (double), (6-34) 

for the optically thin case, and 



i 2u>2p» ± ^a2yi(2tf;2 + 1) - 2{wp^ f 

and A = -/3*+p* (double), (6-35) 
for the optically thick case with = aV/w and /(,^) = 0. 

§7. Hydrodynamics 

A conservative form of the hydrodynamic equations is derived from 

V„(pn") = 0, (7-1) 
7/3.V„(Tfl^1d + rrfd) = 0, (7-2) 
"/3Va(Tfltd + r,fd) = 0. (7-3) 

The first, second, and third equations are the continuity, Euler, and energy equations, 
respectively. For the perfect fiuid, 

T^!^,, = phu'^u^ + Pg'^^, (7-4) 

where h is the specific enthalpy defined by 1 + e + P/ p, and e and P are the specific 
internal energy and pressure. 

For the neutrino-radiation hydrodynamics, it is further necessary to solve the 
continuity equation for leptons or electrons: The continuity equation for electrons is 
written in the form 

VaipYeU'') = pQe, (7-5) 
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where Yg denotes the electron number per nucleon and Qe the electron generation 
rate, determined by the electron capture and neutrino capture by nucleons or nuclei. 
The explicit forms for these equations are 

dtiVlPw) + dji^pwv^) = 0, (7-6) 

= V7 [ - Ph^^a + jkd.p'' + ^S^^dajk - aS-7ai] , (7-7) 

= a^iS'^Kij - f'^jidk In a + 5°n,] , (7-8) 
dt{,/^pwYe) + dji^pwYeV^) = pQeaV^, (7-9) 

where = u'^/u^, and 

ji = -T^^i^nal^i = pwhui, (7-10) 

= T'flu^d"'""'/? = Ph^'^ - -f' (7-11) 
S'" = T^L^Jlf = phV^V^ + P^^K (7-12) 

As in the equations for the radiation moments {E(^^'j, F^^-j^), the Euler and energy 
equations have conservative forms. We note that extension to the magnetohydrody- 
namic equations is straightforward (e.g., Ref. I18p). 

To guarantee the conservation of total momentum and energy, it may be useful 
to solve 

^t[Vlij^ + m + ^,[Vl{j^v' + a{P5,' + P^') - 13^ Fi}] 

= Vl[-iPh + E)dia + (ifc + Fk)di(3^ + + P^'')dajk] , (7-13) 

dt[^{p^ + E)] + d, [^{pW + P{v^ + P') + aF^ - f3^E}] 

= a^[{S'^ + P'^)Ki, - f'^iji + F,)dk In a]. (7-14) 

We note that E, F^, and P^^ here denote the sum of the contribution from all the 
neutrino species. 

§8. Slow-motion limit 

Here, we derive the radiation hydrodynamics equations in the case that (i) the 
spacetime is flat and (ii) the typical velocity of the matter field is much smaller 
than the speed of light. These approximations are often used in the radiation hy- 
drodynamics with Newtonian gravity. Here, several additional words are necessary 
to clarify the condition (ii). First, we denote the typical time and length scales for 
the variation of the matter field by T and L, respectively, and the velocity by V. 
Then, the order of T is equal to L/V, and the acceleration and shear of the matter 
are of order V/T ^ V"^ /L. In the Newtonian approximation for the radiation hy- 
drodynamics, we take into account all the terms associated with the first order in 
V relative to the lowest-order term, but neglect the terms more than second order 
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in V; terms of 0{V'^) such as acceleration and {diVj)"^ are neglected. Because the 
Newtonian potential is the quantity of order V"^, we also neglect the contribution by 
this in the radiation moment equations. 

Then, the equations for the radiation moments defined in the fluid comoving 
frame, (J(^), -H'(^)), are 

d{Lf\diVj) 

dtJi.) + d,iJ^u)v' + H^t)) - V = K(,)(J(^^^) - J(,)), (8-1) 

where we set n° = 1 + O(y^), no = -1 + 0{V'^), = Ui = = Vi, a'' := u^^V ^u^ = 
0{V^), and 



+ ^^V^ (^H^,^^' + ^(i^.^^ + H^^^d.v,) , (8-3) 

For simplicity, we here take into account only the neutrino emission, absorption, and 
isoenergy scattering. The derived forms of Eqs. (jS-ip and (j8-2p agree with those in 
the standard textbooks (e.g., Ref. [B])). On the other hand, the equations for the 
radiation moments defined in a laboratory frame, Fi^^-^j), are 

d{vL,^\diVj) 

dtE^,) + ^£ = ( - + (8-5) 



dv 



where 



+ (^(J)^.^^' + %d,v^ + , (8-7) 

L;{ = 3x-l-^HfH^M + l^i,, .5^., (8.8) 



25 



and we used 



•10) 
•11) 



Again, we note that v is the frequency in the fluid rest frame (not in the laboratory 
frame). As expected, Eqs. (|8-5|) and (|8^6|) have a conservative form. 
The hydrodynamic equations are 



dtp + dj{pv^) = {), 

dt{pvi) + djipviv^ + P6-^) = -pdicj)^ 

1 



+ 



5, 



p[e 



pe + P + -pv^jv^ 



dt{pYe)+dj{pYeV^) = pQe, 



•12) 
•13) 

•14) 
•15) 



where (p-^ is the Newtonian potential. We note that when taking the Newtonian hmit 
of general relativistic hydrodynamics equations, the conservative rest-mass density 
pw^ is replaced to p. The total energy equation is written as 



p(e + y ) +^ 



pe ^ P +\^pv^\v^ + 



-pv'd^(j)N. {8-16) 



§9. Summary 



We derived a truncated moment formalism for general relativistic radiation hy- 
drodynamics modifying the Thome's original formalism.l^ The equations for the 
radiation field are written for the variables defined in the laboratory frame as well 
as in the fluid local rest frame, although the argument angular frequency for the 
radiation moments is always the frequency measured in the fluid local rest frame. 
In the former case, the equations are written in a conservative form (for E'^^) and 
F(iy)i) and essentially the same as those for the hydrodynamic equations in general 
relativity. Thus, they seem to be useful for a well-resolved numerical simulation. 

The source terms are written, focusing on the neutrino transfer in the assump- 
tion that anisotropy of the scattering kernel is small. Then, a formalism for the 
radiation hydrodynamics in numerical relativity is derived in a closed form, assum- 
ing a physically reasonable closure relation among the radiation stress tensor, energy 
density, and energy flux. As long as the radiation fleld is not extremely anisotropic 
in the fluid rest frame (in the optically thick medium) , the employed approximation 
should work well. One merit in the present formalism is that we do not have to 
perform any coordinate transformation when computing the source term, because 
the angular frequency for the radiation fleld is deflned in the fluid local rest frame. 
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The derived equations constitute wave equations for the radiation field. The 
closure relation and variable Eddington factor are appropriately chosen so that the 
characteristic speed is smaller than the speed of light in the free-streaming and grey 
regions. We also notice that (i) for the derivation of the basic equations for the 
radiation field, we do not assume that the fluid velocity is much smaller than the 
speed of light, and (ii) with the chosen closure relation, the effect associated with the 
fluid motion (fluid expansion, acceleration, and shear) may be taken into account. 
Thus, the derived formalism can be employed for the radiation field associated with 
a fast motion, e.g., a fluid moving in the vicinity of a black hole. 

In this formalism, we need to solve 3+1+1 equations (3 is space, 1 is time, and 
frequency space) of F(^)j) or 3+1 equations of {E,Fi) for the radiation part. 

For both cases, the equations for these variables are written in a conservative form, 
and hence, conservation of mass and moments is likely to be well achieved in the 
formalism with these quantities. For the 3+1+1 case, the equations, including ab- 
sorption, emission, and scattering term for neutrinos, are written in the closed form, 
and hence, we do not have to assume anything further. For the 3+1 case, compu- 
tational costs will be saved significantly, but we have to impose several additional 
conditions when performing the frequency integral for the source term: We have to 
assume certain functions for Jj-j^-j and if^"^ (e.g., Ref. ED), for which a physically 
appropriate assumption is required. 

The truncated moment formalism may be a starting point for upgrading the 
current leakage scheme for general relativistic radiation hydrodynamics (e.g., Ref. l20p 
for a review). The leakage scheme is often used for phenomenologically incorporating 
radiation cooling and for a relatively inexpensive radiation hydrodynamic simulation. 
The method is usually quite phenomenological: One first determines optically thick 
and thin regions, respectively, using a rather approximate prescription. Then, for 
the optically thick region, one assumes that the radiation escapes in a diffusive 
manner and for the optically thin region, the radiation escapes freely. In general 
relativistic leakage schemes ,'23 one incorporates the cooling effect in the right-hand 
side of the hydrodynamics equations {VaT'^^ = —S^), and in addition, an equation 
for radiation four-vector field is evolved for the optically thin region. Namely, the 
basic equations are quite similar to those derived in this paper. What is different 
is on the treatment for the source term 5'^^*; in the leakage scheme proposed so 

far, the source term is determined in a quite phenomenological manner. If 5'^^* 
is approximated in a more strict manner starting from the basic equations derived 
in this paper, it will be possible to derive a better and well-funded leakage scheme. 
Furthermore, it will be possible to incorporate the frequency-dependent effect as well 
as neutrino heating. Such work is left for the future. 

Finally, the truncated moment formalism derived here may be used for the trans- 
fer of photons by exchanging the source terms (if we may assume that anisotropy of 
the scattering kernel is small). For example, this formalism will be useful for study- 
ing an accretion flow in the vicinity of stellar-mass and supermassive black holes in 
general relativity. 
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Appendix A 

Stationary radiation in the spherical dilute medium 

In this appendix, we show the solution of a radiation field in the Bondi Rov^^ 
composed of a dilute medium (i.e., optically thin approximation is assumed to work) 
and illustrate that gravitational redshift and Doppler effects are appropriately taken 
into account in the radiation spectrum observed at infinity in our formalism. 

As discussed in § [31 in the optically thin region (where = 0), the radiation 
moments may be written by 

H^^^ = j^^^r, = (A-i) 

where denotes a unit spatial vector for which the spatial component is composed 
only of the radial component. Substitution of these relations into Eq. (|3-22|) with 
= yields 

dJ 

V^{J^,)n-v^efVpu^ = Q, (A-2) 

where /° := n"-|-£° is a null vector. Substitution, in addition, of N^'^^'^ = J (^^^i"- t'^ 
into Eq. (j3-23p yields the same equation as Eq. ()A-2p . 

For a solution of the four velocity, n", we here consider a stationary spherical 
accretion flow (Bondi accretion flow) in the spacetime of a spherical black hole of 
mass M. We choose the line element in the Kerr-Schild coordinates. 



ds^ = -[} ~)'^'^ + —dtdr + (1 + —;—)'^^ + + Od<f'l^-?,) 

in which the coordinate singularity at r = 2GM does not give any messy problem. 
We note however that an analytic solution is also easily derived in the Schwarzschild 
coordinates. 

We denote the infall velocity by u'^ = —u{r) < (cf. for the wind solution 
u''' > 0), and note the relations it = u and £r = n*, where —ut = u*(l — 2GM/r) + 
2GMu/r = y^l -hn2 - 2GM/r. Thus, F = -u + y^l + - 2GM/r > 0, and 
= F{r + 2GM)/{r — 2GM). Using these relations, we flnally reach 

dF 

ri^Vpu^ = ^. (A-4) 

Substitution of this relation into Eq. (|A-2p gives 

^,dJ(^^J^d{r^^^dJ^dr_^ 
dr dr du dr ' 
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and assuming that J^^,) > and F > 0, we obtain 



(A-6) 



dr du dr 

where y(^-^^ := Jiy-^Vr"^. For the case that V is a monotonicahy increasing function of 
r (this is the case for the typical problem), Eq. (|A-6p is rewritten to give 

d\nV d\nv' ^ ' 

Thus, y(^y){r) constitutes a wave equation for the arguments (lnr,lni/), and there- 
fore, the general solution can be derived as 

y(,) = F(/V), or J^^^ = l^, (A.8) 

where F[x) is an arbitrary function of x. Here, — t- 1 for r — )■ oo. Namely, at 
infinity, the observed spectrum is 

J(.) = (A.9) 

On the other hand, for the finite value of r, is smaller than unity. This implies 
that the radiation spectrum should be homogeneously (irrespectively of the value 
of v) shifted to the lower frequency side during outgoing propagation. The redshift 
factor is given by T; an observed radiation with frequency v at infinity is originally 
emitted at a finite radius, remit, with frequency femit = ^/^^{fcmit) > t^- 

Equation (|A-8p indeed captures gravitational redshift and Doppler effects. This 
is clearly found by taking the slow-motion and weak-gravitation approximation for 
F as 

GM 

r«l-n + — . A-10 

2 r 

The first, second, and third terms denote the Doppler, second-order Doppler, and 
gravitational redshift effects, respectively. 
Integration of J(^y^ by v gives 

(irr) 



J= I di^Ji^,) = (A-11) 



where Lq denotes a total flux 

Lo = / duF{u). (A-12) 



Thus Fr may be regarded as a luminosity distance. Note that for r — )• 2GM, V — )• 0. 
Thus, J — )• cxD at r = 2GM; the solution is similar to a solution in the fiat spacetime 
in which a "point" source exists at origin (here which is located at Fr = 0). 

Finally, we point out that the solution derived here will be used as a test-bed 
problem for checking the reliability of a radiation hydrodynamic code based on the 
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truncated moment formalism (note that with Eq. (jA-ip the closure relation holds). 
We also note that the solution given here holds not only for the Bondi flow but also 
any solution in the stationary, spherically symmetric spacetime as long as F is a 
monotonic function of r. 



Appendix B 

Radiation flow in the spherical dilute medium 



Next, we analyze a time-dependent spherically symmetric radiation flow in the 
Schwarzschild spacetime. Again, we assume that neutrinos propagate in the optically 
thin medium. The purpose of this section is to clarify a nature of the closure relations 
(j6-3p and (j6-5p . For this, we ignore the frequency-dependent effects and analyze 
Eqs. ()3-39p and ()3-40p . For the background metric, we again adopt Eq. ()A-3p . In 
this case, the necessary geometric quantities are 

/ 2M\-i/2 2M 2M 

V r / r + 2M r 

2M(r + M) , , 

where we use the units of G = 1 (or we may say that GM is replaced to M). 
Because of the spherical symmetry, we only need to consider the radial component 

1/2 

of radiation moments, F'^ . For the following, we define F := F'^^^r ■ Then the 
condition g^uTj^t^^ = for the closure relations ()6-3p and (j6-5p is written as 

E = \F\ (B-2) 

For the closure relation (j6-3p . the equations for E and F are 

2M , ^ f 2M(2r + M) 3M 

r + 2M'' ^ r + 2M^ ^ r(r + 2M)2 + (r + 2M)2"' ^ 
■ 2M r , 2M{2r + M) 3M 
^ ~ r + 2M^ ^ r + 2M^ ^ r(r + 2M)2 (r + 2M)2^~ ' 

where e = Er'^jrl'^ and / = Fr'^'yrr'^ . The dot (e) and dash (e') denote dte and dre, 
respectively. Defining u± = e ziz f , we obtain two independent equations 

r-2M , M{7r + 2M) 

u+ H u', H 7^ -T^u+ = 0, (B-5) 

+ r + 2M + r(r + 2M)2 + ' ^ ^ 

ii^-u' + , ^ -n_ = 0, (B-6) 

~ r{r + 2M) ' ^ ^ 

and e = (?x+ + U-.)/2 and / = (ti_|_ — n_)/2. This implies that in the absence of u+ 
or u_, the condition ()B-2p is satisfied, but in general, it is not. In particular, for 
the point which satisfies / = {uj^ = u„), one of the characteristic speed becomes 
infinity (see Eq. (j6-29p ). 1*^1 Thus this closure relation should be prohibited for such 

*' In the analysis of spherically symmetric flow here, the characteristic speeds are (r — 2M)/(r + 
2M) and —1. However, if we solve the equation in the Cartesian or cylindrical coordinates, the 
extra characteristic speed (|6-29p appears. 
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a situation (this is resolved in an appropriate choice of the variable Eddington factor 
shown in §6.3). 

For the closure relation (|6-5p . the equations for e and / are 

. _ 2M , r 2M 3Mr + 2M(r + M)s 

r + 2M''^r + 2M^^{r + 2M)^''^ r{r + 2M)^ ' > 

■ 2M-rs 2M{M + 2r + rs) M ^ 

^ r + 2M^^ r{r + 2MY ^^{r + 2M)^ ' ^ ' 

where s = 1 (—1) for F^' > (< 0). Defining u = e — sf = {E — sF)r'^'yll'^ , we 
obtain 

2M , 2M-sM , , 

u ttU + tttttU = 0, (B-9) 

r + 2M (r + 2M)2 ' ^ ^ 

• 2M-rs ^, ^ M{2M + 4r + 3rs) ^ , M _ 

^ r + 2M^^ r{r + 2MY ^^{r + 2MY ^ ' 

Thus, there are also two components: One is determined by u which is a mode of 
physically zero characteristic speed because the coefficient of the transport term is 
equal to — /S*". The other is associated with /, which is an outgoing or ingoing mode 
and obeys the same equation as that of u+ and U- for / > and / < 0, respectively. 
u is regarded as an unphysical mode because it is the measure of deviation from the 
condition (]B-2P : if u = is satisfied, we can follow only the physical mode, but this 
will not be in general the case, in particular in the near zone. The important fact, 
however, is that u does not propagate outward and damps exponentially with time 
in the absence of the source term. This implies that in the zone distant from the 
source, u will be zero because the emission source should be zero there. Thus, in the 
distant optically thin zone, the condition ()B-2p is likely to be satisfied. 

It will be useful to give the solutions of Eqs. ()B-5P and ()B-6P : The general 
solutions for these are written as 



U4 



U- 



^ ' 9+{t-n), (B-ll) 



L (r - 2M) 
r 



-r + 2M 

where g± are arbitrarily functions and r* is a retarded time 



1/2 

g-{t-r), (B-12) 



f r + 2M , /r-2M\ 

n = Jdr—^=r + 4Mln{-^). (B.13) 

Appendix C 

— Numerical experiment for free evolution 



To confirm that the closure relation and variable Eddington factor ()6-25p de- 
scribed in §6 work well in the optically thin medium, we numerically solve radiation 
field equations (|3-39p and (|3-40p on a Bondi flow of a Schwarzschild spacetime. The 
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closure relation is written as 

P'^ = ^ — -E \ , + ^ ^' r ^ + WV^ + WV'\, (C-1) 



where x is assumed to be a function oi F = \F\/E = \J ^^iF^F^ /E. The source 
terms are set to be zero {S" = 0) for simplicity. Numerical simulation was per- 
formed assuming the axial and equatorial plane symmetries. As in Appendix A 
and B, the Kerr-Schild coordinates are adopted, and the same Bondi solution as in 
Ref. [T8|) is employed. The basic equations are essentially the same as those solved in 
general relativistic hydrodynamic simulation. We employ the same scheme as used 
in Ref. [T8|) for a solution of E and F^. Specifically, the transport term is handled 
using a Kurganov-Tadmor schemd^D with a piecewise parabolic reconstruction for 
the quantities of cell interfaces. The fourth-order Runge-Kutta method is employed 
for the time integration. The characteristic speed is not analytically computed for 
the general form of Pij. Thus, we simply write it in the linear combination form 

— ^ ^thin H ^ -^thick- V^''^) 

For the case that the relation, E < \F\, is accidentally realized at a point, we set 
X = 1, and Athin is limited to be smaller than unity. 

With the time evolution, the radiation fields flow away from the computational 
domain. To handle this correctly, an outgoing boundary condition is imposed for the 
outer boundaries, and inside the radius r < l.SAf, we artificially set E = Fi^. = 0. 

First, we consider the solution for E = \F\ derived in Appendix B. In this case, 
F is always unity, and thus, x = 1 always holds. For the outgoing flow F = E, the 
solution is written as 



E 1/2 ^ J, 1/2 ^ 1 
^ irr ^ irr ^ 



and for the ingoing flow F = —E, 



(r + 2Mf 
r3(r - 2My 



^1/2 

9+{t-r.), (C-3) 



E^U^ = -F^W = lir'ir + 2M)]-^/^g.{t - r). (C-4) 
We choose a form of a wave packet as 

g+{r,) = exp[-(r, - u^f /m\ g.{r) = exp[-(r - r^f /SM\ (C-5) 

Numerical simulations were performed for r*o = r^{r = 6M) and tq = 32M. The 
computational domain covers a region [0 : 40M] both for x and z with a uniform 
grid spacing O.IM. 

Figure 1 plots the evolution of Ey^i'^{r — 2M)^(1 -|- 2M/r)^^/^ for the outgo- 
ing solution (left) and £^7r/^r^/^(r -|- 2MYI'^ for the ingoing solution (right). This 
shows that besides a small phase error, the numerical solutions reproduce the exact 
solutions. 
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5 10 15 20 25 30 35 40 5 10 15 20 25 30 35 40 

x/M x/M 



Fig. 1. Evolution of outgoing (left) and ingoing (right) solutions. For the outgoing solution, the 
profiles of E-yW^ir - 2M)^(1 + 2M/r)-^/'^ at t = and t = 35M along the x and z axes are 
plotted. For the ingoing solution, the profiles of E-fW^r^^'^{r + 2M)^''^ at t = and t = 25M 
along the x and z axes are plotted. The solid and dashed curves denote the numerical and exact 
solutions. The numerical solutions for x and z axes agree approximately. 



We also performed a simulation for F = at t = with 

^7^1/2 ^ exp[-(r - 10Mf/8M^]. (C-6) 

In this case, x = 1/3 at t = 0, but with the time evolution, \F\ becomes nonzero 
and X becomes larger than 1/3. Because the variable Eddington factor x is varied 
with the evolution, we do not have the exact solution. The purpose is to test if our 
formalism allows a stable numerical solution. 

Figure 2(a) plots the time evolution of the wave packet. Here we plot Ejrr'^r'^ 
along X and z axes (two results approximately agree and cannot be distinguished 
in the figure). After the evolution starts, the wave packet is split into outgoing 
and ingoing parts. Both modes propagate smoothly with no trouble. Figure 2(b) 
plots the evolution of F/E. This is initially zero, but with the free propagation, 
it approaches to unity: At t/M = 5, 20, and 50, F/E is larger than 0.8 for 15 ^ 
x/M <j 25, 15 ^ x/M ^ /35, and x/M ^ 5, respectively. No problem is found for 
the propagation, and thus, as far as the numerical issues are concerned, the closure 
relation and variable Eddington factor employed here have no problem. 
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